Many-body GW Calculations of Ground-State Properties: Quasi-2D Electron Systems 

and van der Waals Forces 
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We present GW many-body results for ground-state properties of two simple but very distinct 
families of inhomogenous systems in which traditional implementations of density-functional theory 
(DFT) fail drastically. The GW approach gives notably better results than the well-known random- 
phase approximation, at a similar computational cost. These results establish GW as an superior 
alternative to standard DFT schemes without the expensive numerical effort required by quantum 
Monte Carlo simulations. 



71.15.Nc,71.15.Mb,71.45.Gm 

Ab initio many-body-theory (MBT) methods, particu- 
larly those based on Hedin's approximation have 
been used extensively during the last decade to calculate 
excited-state properties of electron systems in solid-state 
physics 0]. In addition, there is a recent and increas- 
ing interest in the application of such methods to ob- 
tain ground-state properties f^Tf!- In MBT, electron- 
electron correlations are taken into account directly with- 
out resorting to the mean-field density-based approxima- 
tions used in routine implementations of Density Func- 
tional Theory (DFT) Q, thus providing a more micro- 
scopic description of the interacting many-electron prob- 
lem. 

Modern MBT calculations use efficient algorithms for 
the evaluation of MBT quantities [such as the dielec- 
tric function e(w), the self-energy operator S(w), and 
the one-particle Green function G (u>)] required in the 
study of excited-state properties of real materials JTHpH . 
The application of these techniques to the calculation of 
ground-state properties, such as the density or the to- 
tal energy, opens the appealing possibility of treating all 
the electronic properties in the same fashion. Although 
more expensive than DFT, these MBT methods do not 
require the demanding computational effort of quantum 
Monte Carlo simulations OS. Nonetheless, ground-state 
calculations based on MBT must be painstakingly as- 
sessed. First, approximations that have proved success- 
ful for spectral properties might not necessarily be good 
methods for structural properties. Second, the energies 
and lifetimes of quasiparticles are mainly determined by 
the pole structure of G (w) , while ground-state proper- 
ties emerge from a multidimensional integration of G (ui). 
As a consequence, features that can be safely ignored 
in the determination of quasiparticle properties, like the 
high-frequency behavior or high-energy-transfer matrix 
elements, play an important role if we want to calculate, 
for instance, the ground-state energy. The development 
of general-purpose procedures requires a careful study of 
the optimal treatment of these points. 



In this Letter we present nonselfconsistent (in the sense 
described below) GW calculations for the ground-state 
properties of two families of simple inhomogeneous sys- 
tems: the quasi-two-dimensional (2D) electron gas and a 
pair of interacting jellium slabs. These provide a suitable 
test cases for an extensive assessment of the performance 
of GW-MBT ground-state calculations. First, despite 
the simplicity of these systems, LDA and GGA fail to 
describe them properly. For the quasi-2D gas, the high 
inhomogeneity of the density profile along the confining 
direction is clearly beyond the scope of local or semilo- 
cal approaches fl5|] . The situation is quite different if we 
consider interacting jellium slabs. In the limit of large 
separation, in which the densities of each subsystem do 
not overlap, dispersion or van der Waals (vdW) forces are 
evident. These forces are due to long-ranged Coulomb 
correlations and, hence, cannot be described at all by 
the LDA or GGA HOI. As shown below, even non- 
local DFT prescriptions for the XC energy, such as the 
weighted density approximation (WDA) p7| , are unable 
to reproduce such forces. On the other hand, these sys- 
tems exhibit translational invariance along the xy plane, 
thus reducing significantly the number of independent 
variables needed to describe the spatial dependence of 
all the operators, assisting the production of the highly- 
converged test calculations required here. 

In Hedin's GW framework, the self-energy E of a sys- 
tem of N electrons under an external potential v ex t (r) is 
approximated by 



E(l,2)=tG(l,2+) W(l,2), 



(1) 



where the labels 1 and 2 symbolize space-time coordi- 
nates. W is the screened Coulomb potential, which is 
exactly related to the bare Coulomb potential w and 
the polarizability P by W (u>) — w + wP (u>) W (u>) 
< (lS) w (the usual matrix operations are implied). Un- 
der the GW approach, the polarizability is given by 
P(l,2) = -2iG (1,2) (7(2,1+). Finally, E and G are 
linked through the Dyson equation G _1 (cj) = uj — [t + 



1 



V e l + E (w)], with i the one-electron kinetic energy, and 
W e i ( r ) the classical electrostatic potential [the sum of 
^ext (r) and the Hartree interaction potential]. This set 
of equations defines a selfconsistent problem, but routine 
GW calculations concerned with the quasiparticle prop- 
erties do not attempt selfconsistency: when evaluating £, 
the interacting Green function G is generally substituted 
by that corresponding to the noninteracting Kohn-Sham 
(KS) system under the LDA (or GGA) approximation. 

It is well known that full selfconsistency implies a wors- 
ening of the description of spectral properties p8| . How- 
ever, for the 3D |Q and 2D Q homo geneous elec- 
tron gases, the total energy arising from G by using the 
Galitskii-Migdal formula of MBT |§ fits the (essentially 
exact) QMC values extremely well if G is obtained self- 
consistently. Nonetheless, at the nonselfconsistent level 
(which we will call Go Wo) the exact exchange energy 
is already built in and, at the same time, correlation ef- 
fects beyond the random-phase approximation (RPA) are 
taken into account, as reflected in good total energies and 
very good total energy differences |^Jio|| for these homo- 
geneous systems. Since the relevant quantity for struc- 
tural properties is the total energy difference, Go Wo (and 
even, in some cases, RPA ^0|) is likely to fulfill a use- 
ful role. The known particle-number violation under the 



Go Wo approximation 21 1, while non-zero, is so small that 
it can be ignored in practical applications p^ |. 

All the calculations reported in this Letter are done 
as follows. First, we perform a standard LDA-KS cal- 
culation, in which, owing to the translational invariance 
along the xy plane, the KS orbitals are organized into 
subbands. That is, each KS state is given by 0„k (r) = 



ip n (z) e lk ' p /-\/27r, where k = (k Xl k y ) and p = (x,y) 
denote the two-dimensional momentum and position in 
the xy-plane. Then, the noninteracting polarizability 
Pa (1, 2) = -2iG (1, 2) G (2, 1+) is calculated at imagi- 
nary frequencies using the expression 

N occ oo 

P (z, z', k; iu>) = ^ X! Snmk 7 ™ m ( z ) lnm ( 2 ) 



Here, N occ is the number of occupied KS subbands, 
lnm(z) = ip n (z)ip m (z), and the coefficients S nmk (iuj) 
admit a fully analytical expression pj23|| . Pq(\uj) could 
be also calculated fully numerically [|l 3f| but due to the 
high symmetry of the system, the present method is more 
efficient. The infinite sum in Eq. || is truncated at a cer- 
tain value -ZVb, which acts as a convergence parameter. 

The next step is the inversion of the RPA dielectric 
function eo (iw) = 1 — wPo (iuj) and the evaluation of 
Wo (ioj) = e^ 1 ^)"!!). This can be done in a double 
cosine basis representation but the high inhomo- 

geneity of the systems along the z direction suggests a 
different method based on Ref . |l^] . Note that the polar- 
izability spans a Hilbert space made up by the product 
states 7nm- Obviously, the set {-jnm} has a large num- 



ber of linear dependences which can be eliminated us- 
ing a standard Gramm-Schmidt procedure, so obtaining 
an optimized basis set B op = {( a (z)}. In all the cases 



here studied, N n 



N B + N 



1 ^-functions ensure 



an almost perfect representation (with a relative error 
less than 10~ 4 ) of all the product states 7„ m . Hence, 
P (z, z', k; iui) = E Q ,/3Ca {z)Pq i k ,^) a0 Cfi {*'), where 
the matrix elements Pq (fc, iu) a p can be immediately ob- 
tained from S nm k and the scalar products (jnmXa)- 
Then, we calculate the representation w (k) a/3 of the bare 
Coulomb potential, and the matrix elements W (fc, i^) Q ^ 
of the screened Coulomb potential are easily obtained 
by a matrix inversion for each value of k. Well con- 
verged results are typically obtained with N& — 80 sub- 
bands, although for quasi-2D systems Nb is significantly 
less. RPA correlation energies, where required, can be 
straightforwardly obtained in this representation. 

The real-space representation of Wo is given by the 
expansion 



W (z,z',;ir) =i £ / 



d0J ; 



i(ujT+k-p) 



2tt 



Xt a (z)(p(z')W(k,kj)a0 > 

where ^ k = (27r) -2 J dk, whereas the KS Green function 
Gq [z, z', p; it) is calculated directly from the KS eigen- 
states Q . By using ([!]) we readily get the self-energy op- 
erator in real space and imaginary time and, eventually, 
its representation £ (fc, ir) nm in the KS basis set. We use 
typically 100 (or fewer) t/ui points in a Gauss-Legendre 
grid with ui max — 80A C (A c is the width of the conduc- 
tion band), and 100-150 p/k points with fc^ ax /2 ~ a7 max . 
It is important to emphasize that the asymptotic time 
and frequency tails of all the operators must be treated 
analytically to ensure smooth and rapid convergence. 

In the imaginary-time/frequency representation 
£ (/i + \u)) = — i J dr E (ir) e~ la,T , with p the interacting 
chemical potential. Hence, the Green function at imagi- 
nary frequencies is the solution of the Dyson equation 



G 1 (p + iuj) = p + iui — h 



S( M . 



IUJ 



vxc 



(3) 



with /iks — t + vs the KS-LDA hamiltonian and vxc (r) 
the LDA-XC potential. G is calculated in the KS rep- 
resentation G (k,iu) nml with p previously obtained by 
diagonalizing the quasiparticle hamiltonian at io; = 0, 
^QP (m) = ^ks + £ (p) — vxc, an d by imposing that the 
volume enclosed by the interacting Fermi surface equals 
the KS value. A small term Sv(r), accounting for the 
change in the Hartree potential due to the differences be- 
tween the Go Wo and LDA densities, should be included 
into (^) and /iqp- However, this term induces an imper- 
ceptible change in the shape of G, and it has a negligible 
influence on the ground-state properties. 

The electron density and energy are given, respec- 
tively, by n(r) = tt^ 1 § duo G (r, r; iu) and E = E c \ + 
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FIG. 1. XC energy per particle for thin confined jellium 
slabs of fixed 2D density n 2D = 3/4 7r as a function of their 
thickness L. Lines: LDA, WDA, and RPA; circles: GoWo- 
The exchange energy (dash-dotted line) has been included as 
a reference. Note the effect due to the filling of a second 
subband at L ~ 4.2ao, which is especially evident in ex- The 
arrow indicates the QMC value for the 2D limit. 

7T _1 tv § duj[i + |S (ilo)]G (iuj), where E c \ is the classi- 
cal electrostatic energy [which can be obtained from 
n(r)], and tr symbols the spatial trace. Owing to the 
symmetry of our model systems (the generalization to 
arbitrary geometries is straightforward), we can write 

n (z) = n a (z) + Sn (z) = n (z) + J2 n ,k iV'n (z)| Sf nk , 
no (z) being the LDA density, and 

/d 
— [G (k, n + iu) nn - G (k, mo + iw) nn ] , (4) 

(fio is the LDA chemical potential). On the other 
hand, the exchange energy per surface unit is given 
directly by E x /S = \ J2 n ,k /nfe^x where E x 

is the frequency-independent part of the self-energy 
£x(ri,r 2 ) = ilim r ^ + G (ri,r 2 ;iT)w(ri 2 ) (i.e., the 
Fock operator). Finally, the correlation energy Ec is 
written as the sum of its kinetic-energy and electrostatic 
parts: 



= 51 SnkSfnk ~ J dzv s (z) 5n {z 

n,k 



(5) 



Wc = 



n ,k 



2 E x (k) nn 5f nk + 



/duj 
— £ c (fc, (i + iuj) nm G (fc, n + iw) r 



(6) 



where we have included the frequency-dependent part of 
the self-energy £c (iw) = £ (iu>) — Ex- As usual, we sup- 
pose that the total electrostatic energy is given correctly 
by the LDA, since we have checked carefully that the 
density variation Sn (z) causes only minor changes in the 
electronic energy. We also note that when evaluating Ec 
we do not need all the matrix elements of £ and G. In 



general, full convergence is achieved by performing the 
sum only over LDA states such that £„k J$ 6A C , and this 
also applies to the resolution of Dyson's equation. There- 
fore, the most demanding part of the calculation is the 
evaluation of e _1 (iu>) which, in any case, is required to 
obtain the RPA correlation energy. 

To mimic a quasi-2D electron system, we have taken 
a thin jellium slab with a background density n = 
(l 777 "^) -1 and width L. This slab is bounded by two 
infinite planar walls, and overall charge neutrality is as- 
sumed. We keep the number of particles per unit surface 
area n 2 D = nL constant, in such a way that the limit 
L — > corresponds to a 2D homogeneous electron gas 
(HEG) with density n 2 D- In Fig. 1 wc plot, for several 
values of L, the XC energy per particle exc given by 
the LDA, the nonlocal WDA, the RPA, and the G Q W 
method. As commented previously, the LDA diverges 
when approaching the 2D limit, whereas the WDA be- 
haves reasonably well, slightly underestimating the ab- 
solute value of exc in the strict 2D limit. The RPA 
does not show any pathological behaviour, but it overes- 
timates | exc | by more than 20 mHa/e~ for all configura- 
tions. Finally, Go Wo, whose superiority to the RPA has 
already been established in the 3D limit, ||[l0) retains 
this in the transition to the 2D limit. Its performance is 
similar to that of the WDA for these systems (although 
the residual error for the 2D gas has the opposite sign). 
It is important to point out that the RPA and Go Wo XC 
energies obtained by starting from a WDA-KS calcula- 
tion are indistinguishable from those plotted in Fig. 1. 
Thus, the specific choice of the KS method seems to be 
of minor importance when calculating XC energies using 
MBT. 

The study of the interacting energy between two un- 
confined jellium slabs is of more direct significance, as 
it has been considered as a benchmark for seamless cor- 
relation functionals attempting to describe vdW forces 
^,|4j. By varying the distance d between the slabs we 
cover configurations in which the density profiles of each 
subsystem overlap (i.e. there is a covalent bond), and 
situations (when d ^> 0) in which there is no such over- 
lap and the only source of bonding is the appearance of 
vdW forces. In the upper panel of Fig. 2 we plot the 
XC energy per particle exc (d) a s a function of d using 
the LDA (the WDA gives very similar results), the RPA, 
and the Go Wo, for two slabs of width L = 12ao and a 
background density n corresponding to r s = 3.93. We 
can see that the GoW) reduces the RPA error by 60%. 
In the lower panel we represent the correlation binding 
energy per particle [defined as ec (d) = ec (d) — ec (oo)]. 
First, neither the LDA nor even the nonlocal WDA is able 
to reproduce the characteristic asymptotic dr 2 vdW be- 
havior. This behavior is described by the RPA and the 
Go Wo, and for both approximations ec (d 3> 0) is very 
similar (which is not a surprise because such an asymp- 
totic behavior is fully described at the RPA level H ) . For 
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FIG. 2. Upper panel: XC energy per particle for two jel- 
lium slabs as a function of the distance d. Lines: LDA, WDA 
and RPA; circles: Go Wo; squares: Go Wo plus LDA- like cor- 
rection term. Lower panel: correlation binding energy per 
particle. The LDA and WDA contributions have been ob- 
tained by subtracting the exact exchange energy from the 
corresponding XC values. The exchange binding energy 
(dashed-dotted line) has been also included in this panel. 

intermediate and small separations, there are slight dif- 
ferences between the RPA and the Go Wo but much less 
important than those appearing when comparing the to- 
tal correlation energies. 

A detail that is worth pointing out is the fact that 
the error in the absolute Go Wo correlation energy is 
amenable to a LDA-like correction. From the differ- 
ences, given in Ref. between QMC and Go Wo cor- 
relation energies for the HEG, we can build a functional 
AE C = J drn(r) Sec ["(r)]. The absolute G W energy 
corrected in this way is in broad correspondence with the 
LDA energy, while the binding energy retains its correct 
1/d 2 behavior at large d and is little altered for small d 
(see Fig. 2). 

In summary, we have performed GW-MBT calcula- 
tions to evaluate the ground-state energy of inhomoge- 
neous systems. With practically the same cost we have 
obtained correlation energies beyond the random phase 
approximation. The importance of a truly ab initio treat- 
ment of electron many-body effects is evident in the 
model systems we have chosen, for which traditional im- 
plementations of DFT are completely inaccurate. 
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